---
title: "Interval Truth Model"
subtitle: "Visualizations of Main Simulation Study"
author:
- name: Matthias Kloft
orcid: 0000-0003-1845-6957
affiliations: University of Marburg
- name: Björn S. Siepe
orcid: 0000-0002-9558-4648
affiliations: University of Marburg
- name: Daniel W. Heck
orcid: 0000-0002-6302-9252
affiliations: University of Marburg
date: "`r Sys.Date()`"
format:
html:
toc: true
number-sections: true
theme: cosmo
code-fold: true
code-tools: true
code-summary: "Show the code"
fig-width: 7
fig-height: 4.5
embed-resources: true
execute:
message: false
warning: false
---
# Background
This file contains the visualizations of the simulation study for the Interval Truth Model.
The results of the preliminary simulation study on the two different link functions can be found in the file `src/01_simulation_study_link_functions_visualizations.Qmd`.
All errorbars in this document represent $\pm 1$ Monte Carlo Standard Errors.
# Prep
We first load all relevant packages:
```{r load-pkgs}
packages <- c(
"tidyverse",
"SimDesign",
"here",
"psych",
"ggh4x",
"ggokabeito",
"ggExtra",
"showtext",
"ggdist",
"pander",
"sysfonts",
"DT"
)
if (!require("pacman")) install.packages("pacman")
pacman::p_load(packages, update = F, character.only = T)
# if(!require("cmdstanr")){
# install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
# library(cmdstanr)
# }
# default chunk options
knitr::opts_chunk$set(
fig.height = 7,
fig.width = 10,
include = TRUE,
message = FALSE,
warning = FALSE
)
source(here("src", "00_functions.R"))
# specify path to full simulation results
path_full_results <- here::here(
"sim_results",
"sim_results_itm_2025-04-28_03-30-27",
"full_sim_results_itm"
)
# add google font
sysfonts::font_add_google("News Cycle", "news")
# use showtext
showtext::showtext_auto()
```
As the consensus locations and width have slightly different standard deviations, we can use these to standardize performance measures for ease of interpretability.
Obtain the SDs from the data-generation function:
```{r}
link <- "ilr"
mean_benchmark <- simplex_to_bvn(c(.4, .2, .4), type = link)
sd_benchmark_loc <- simplex_to_bvn(c(.98, .01, .01), type = link)
sd_benchmark_wid <- simplex_to_bvn(c(.495, .01, .495), type = link)
# mean for Tr_loc
mu_Tr_loc <- mean_benchmark[1]
# mean for Tr_wid
mu_Tr_wid <- mean_benchmark[2]
# SD for Tr_loc
sigma_Tr_loc <- sd_benchmark_loc[1] / 4
# SD for Tr_wid
sigma_Tr_wid <- abs(sd_benchmark_wid[2] - mean_benchmark[2]) / 4
# SD for the mean of Tr_loc and Tr_wid
sigma_Tr_interval <- sqrt(sigma_Tr_loc^2 + sigma_Tr_wid^2) / 2
# SD for a_loc on the log-scale
sigma_a_loc <- .3
# SD for b_loc
sigma_b_loc <- sigma_Tr_loc / 3
# SD for b_wid
sigma_b_wid <- sigma_Tr_wid / 3
# SD for lambda_loc and E_loc on the log-scale
sigma_lambda_E_loc <- .3
# SD for lambda_wid and E_wid on the log-scale
sigma_lambda_E_wid <- .3
```
# Deviations from Preregistration
| Deviation | Explanation |
|-----------|-------------|
| We used `rstan` instead of `cmdstanr`. | We used `rstan` due to its ease of installation and use for applied users, and its useful helper functions. This should not influence the quality of our results. |
| We did not use multivariate priors for proficiency and discernibility parameters. | TODO give reasoning. |
| We set $\omega_j = 0$ in data generation. | We simply forgot to mention this in the preregistration. TODO @Matze, is this correct? |
| Changed formula for absolute bias and MSE, with denominator changing from number of participants to the number of items. | |
| We used 500 instead of 1,000 simulation repetitions | We underestimated the computational workload, but checked the uncertainty of point estimates and realized that it was so small that 500 repetitions sufficed, which is in line with our preregistered criterion of a certain maximum MCSE. |
| We changed the R version. | During the revision, we used the newest R version. This should not change the results.|
| We added additional benchmarks (i.e., simple medians). | We were asked to do so during the revision, and added it because it provides an even stronger test of our model. |
| We added the Pearson correlation between true and estimated parameter values as an additional performance measure. | We were asked to include the performance/recovery of the secondary parameter besides the consensus intervals. The correlation facilitates the interpretation of results regarding recovery. |
Additionally, as we already explained in the preregistration, we did not preregister the full model setup including the priors, which left us with some freedom to choose the priors. We make these choices transparent in the manuscript.
***
# Results
The full results of the simulation study are available in two data frames:
- `sim_res_itm_server.rds`: Results of the simulation study, missing some MSE results due to a small coding error, but containing all information about seed, RAM usage, etc.
- `sim_res_itm_0808.rds`: Resummarized results of the simulation study, now also containing all MSE results as well as more information on divergent transitions.
```{r}
sim_res_itm <- readRDS(
here::here(
"sim_results",
"sim_results_itm_2025-04-28_03-30-27",
"resummarized_results_itm.rds"
)
)
```
Prepare data and convert to long format for plotting:
```{r}
sim_res_itm_long <- sim_res_itm |>
select(!c(contains("conv"), contains("divtrans"), contains("_sd"))) |>
# the following columns aren't needed any more for the "resummarized" results
# "REPLICATIONS", "SIM_TIME", "RAM_USED", "SEED", "COMPLETED", "WARNINGS")) |>
# delete "_fn_" from every column name
rename_all( ~ str_remove(., "_fn")) |>
pivot_longer(
cols = !c(n_respondents, n_items),
names_to = "measure",
values_to = "value"
) |>
# only look at non-log transformed values
filter(!grepl("log", measure)) |>
# rename for easier separation
mutate(measure = str_replace(measure, "abs_bias", "absbias")) |>
# add string so that each measure has same number of underscores
mutate(measure = str_replace(measure, "omega", "omega_cor")) |>
# remove only the first underscore
mutate(measure = sub("_", "", measure, fixed = TRUE)) |>
separate_wider_delim(measure,
names = c("measure", "summary", "pm", "param"),
delim = "_") |>
group_by(n_respondents, n_items, measure, summary, pm) |>
pivot_wider(names_from = "param", values_from = "value") |>
ungroup() |>
mutate(n_respondents = factor(n_respondents))
```
## Convergence and Rhat
We now present the average convergence statistics. All values are averages over all replications. A performance measure with a $>$ or a $<$ indicates how many out of 1000 replications had a certain property. For example, $\hat{R} < 1.1$ means that on average, a certain amount of replications had an $\hat{R}$ value below 1.1.
```{r}
# Renaming for the table
name_mapping <- c(
"rhat_mean" = "$\\hat{R}_{\\text{mean}}$",
"rhat_prop1c1" = "$\\hat{R} < 1.1$",
"rhat_prop1c05" = "$\\hat{R} < 1.05$",
"rhat_prop1c01" = "$\\hat{R} < 1.01$",
"ess_bulk_mean" = "$\\text{ESS}_{\\text{bulk, mean}}$",
"ess_bulk_prop100" = "$\\text{ESS}_{\\text{bulk > 100}}$",
"ess_bulk_prop200" = "$\\text{ESS}_{\\text{bulk > 200}}$",
"ess_bulk_prop300" = "$\\text{ESS}_{\\text{bulk > 300}}$",
"ess_bulk_prop400" = "$\\text{ESS}_{\\text{bulk > 400}}$",
"ess_tail_mean" = "$\\text{ESS}_{\\text{tail, mean}}$",
"ess_tail_prop100" = "$\\text{ESS}_{\\text{tail > 100}}$",
"ess_tail_prop200" = "$\\text{ESS}_{\\text{tail > 200}}$",
"ess_tail_prop300" = "$\\text{ESS}_{\\text{tail > 300}}$",
"ess_tail_prop400" = "$\\text{ESS}_{\\text{tail > 400}}$"
)
sim_res_itm |>
select(contains("conv")) |>
summarize(across(everything(), mean)) |>
# remove "mean_conv." from every column name
rename_all(~str_remove(., "mean_conv.")) |>
rename_all(~str_remove(., "mcse_conv.")) |>
pivot_longer(cols = everything()) |>
# separate mean and mcse based on last underscore
separate(name, into = c("name", "suffix"), sep = "_(?=[^_]+$)", remove = FALSE) |>
pivot_wider(names_from = suffix, values_from = value) |>
mutate(mean = round(mean, 4),
mcse = round(mcse, 4)) |>
mutate(name = name_mapping[name]) |>
knitr::kable()
```
The divergent transitions are shown in the table below. For purposes of readability, we omit the MCSEs here.
```{r}
sim_res_itm |>
select(n_items, n_respondents, contains("divtrans")) |>
select(-contains("mcse")) |>
rename(
"Mean DivTrans (all)" = "mean_divtrans.divergent_transitions_mean",
"Proportion with DivTrans" = "nonz_divtrans.divergent_transitions_nonzero",
"Mean DivTrans (models with any DivTrans)" = "nonzmean_divtrans.divergent_transitionsnonzero_mean"
) |>
mutate(
n_items = factor(n_items, ordered = TRUE),
n_respondents = factor(n_respondents, ordered = TRUE)
) |>
arrange(n_items, n_respondents) |>
rename(Items = n_items, Respondents = n_respondents) |>
# round numeric columns
mutate(across(where(is.numeric), ~ round(.x, 3))) |>
knitr::kable(align = c("c", "c", "r", "r", "r"))
```
Check if our prespecified MCSE criteria ($<.05$) was fulfilled in all conditions:
```{r}
sim_res_itm |>
select(all_of(contains("mcse"))) |>
select(all_of(contains("bias"))) |>
summarize(across(everything(), max)) |>
pivot_longer(cols = everything()) |>
arrange(desc(value)) |>
knitr::kable()
```
## Visualization: Consensus Location + Width
### Absolute Bias of Consensus Location + Width
#### Combined AbsBias
Here, we show the absolute bias of the consensus location and width for the different sample sizes and number of items:
```{r}
plot_sim_absbias <- sim_res_itm_long |>
filter(pm == "absbias") |>
filter(summary == "mean") |>
filter(measure %in% c("Trinterval", "simplemeaninterval", "simplemedinterval")) |>
mutate(
measure = case_when(
measure == "Trinterval" ~ "ITM",
measure == "simplemeaninterval" ~ "Simple Means",
measure == "simplemedinterval" ~ "Simple Medians"
)
) |>
mutate(n_items = paste0(n_items, " Items")) |>
# order n_items correctly
mutate(n_items = factor(n_items, levels = c("5 Items", "10 Items", "20 Items", "40 Items"))) |>
# # Standardize with true variability
# mutate(mean = mean / sigma_Tr_interval,
# mcse = mcse / sigma_Tr_interval) |>
ggplot(aes(
x = n_respondents,
y = mean,
color = measure,
group = measure
)) +
# add vertical line between different sample sizes
geom_vline(colour = "#F3F4F5", xintercept = seq(1.5, 4, 1)) +
geom_point(position = position_dodge(0.7), size = 2.5) +
geom_errorbar(
aes(ymin = mean - 1 * mcse, ymax = mean + 1 * mcse),
width = 0.95,
position = position_dodge(0.7),
show.legend = FALSE
) +
ggh4x::facet_wrap2(n_items ~ ., axes = "all", nrow = 1) +
scale_y_continuous(limits = c(0, .49), expand = c(0, 0)) +
ggokabeito::scale_color_okabe_ito(order = c(5, 1, 3)) +
labs(x = "Number of Respondents", y = "Absolute Bias", color = "") +
theme_icm() +
theme(legend.position = "top", text = element_text(size = 22))
# For manuscript plot:
# theme(legend.position = "top",
# text = element_text(size = 14),
# axis.text.x = element_text(size = 13, margin = margin(b = 10)),
# axis.text.y = element_text(size = 13, margin = margin(l = 10)),
# legend.text = element_text(size = 15.5))
#
ggsave(
here("plots", "sim_main", "sim_absbias.pdf"),
plot_sim_absbias,
width = 9,
height = 4.5
)
plot_sim_absbias
```
#### Separate AbsBias
Here, we focus on the ITM and separate the bias by location and width, while standardizing for true variability:
```{r}
plot_sim_bias_widloc <- sim_res_itm_long |>
filter(pm == "absbias") |>
filter(summary == "mean") |>
filter(measure %in% c("Trloc", "Trwid")) |>
mutate(measure = case_when(
measure == "Trloc" ~ "Location",
measure == "Trwid" ~ "Width"
)) |>
mutate(n_items = paste0(n_items, " Items")) |>
# order n_items correctly
mutate(n_items = factor(n_items, levels = c("5 Items",
"10 Items",
"20 Items",
"40 Items"))) |>
# Standardize with true variability
mutate(mean = if_else(measure == "Location",
mean / sigma_Tr_loc,
mean / sigma_Tr_wid
),
mcse = if_else(measure == "Location",
mcse / sigma_Tr_loc,
mcse / sigma_Tr_wid
)) |>
ggplot(aes(x = n_respondents,
y = mean,
color = measure,
group = measure)) +
# add vertical line between different sample sizes
geom_vline(colour = "#F3F4F5", xintercept = seq(1.5, 4, 1))+
geom_point(position = position_dodge(0.7),
size = 2.5) +
geom_errorbar(aes(ymin = mean - 1*mcse,
ymax = mean + 1*mcse),
width = .8,
position = position_dodge(0.7),
show.legend = FALSE)+
ggh4x::facet_wrap2(n_items ~ .,
axes = "all",
nrow = 1) +
scale_y_continuous(limits = c(0, 0.31), expand = c(0,0)) +
ggokabeito::scale_color_okabe_ito(order = c(2, 7))+
labs(x = "Number of Respondents",
y = "Standardized Absolute Bias",
color = "") +
theme_icm()+
theme(legend.position = "top",
text = element_text(size = 28))
ggsave(here("plots","sim_main", "sim_absbias_widloc.pdf"),
plot_sim_bias_widloc, width = 7, height = 3.75)
plot_sim_bias_widloc
```
Alternatively, combine plots for location and width and standardize with true variability:
```{r}
plot_sim_absbias_widloc_comb <- sim_res_itm_long |>
filter(pm == "absbias") |>
filter(summary == "mean") |>
filter(
measure %in% c(
"Trloc",
"Trwid",
"simplemeanloc",
"simplemeanwid",
"simplemedloc",
"simplemedwid"
)
) |>
# split based on model
mutate(
model = case_when(
grepl("^Tr", measure) ~ "ITM",
grepl("^simplemean", measure) ~ "Simple Means",
grepl("^simplemed", measure) ~ "Simple Medians"
),
measure = sub("^Tr|simplemean|simplemed", "", measure)
) |>
mutate(measure = case_when(measure == "loc" ~ "Location", measure == "wid" ~ "Width")) |>
# Standardize with true variability
mutate(
mean = if_else(measure == "Location", mean / sigma_Tr_loc, mean / sigma_Tr_wid),
mcse = if_else(measure == "Location", mcse / sigma_Tr_loc, mcse / sigma_Tr_wid)
) |>
mutate(n_items = paste0(n_items, " Items")) |>
# order n_items correctly
mutate(n_items = factor(n_items, levels = c("5 Items", "10 Items", "20 Items", "40 Items"))) |>
ggplot(aes(
x = n_respondents,
y = mean,
color = model,
group = model
)) +
# add vertical line between different sample sizes
geom_vline(colour = "#F3F4F5", xintercept = seq(1.5, 4, 1)) +
geom_point(position = position_dodge(0.7), size = 2.5) +
geom_errorbar(
aes(ymin = mean - 1 * mcse, ymax = mean + 1 * mcse),
width = .8,
position = position_dodge(0.7),
show.legend = FALSE
) +
ggh4x::facet_grid2(measure ~ n_items, axes = "all") +
scale_y_continuous(limits = c(0, 0.39), expand = c(0, 0)) +
ggokabeito::scale_color_okabe_ito(order = c(5, 1, 3)) +
labs(x = "Number of Respondents", y = "Standardized Absolute Bias", color = "") +
theme_icm() +
# theme(legend.position = "top", text = element_text(size = 14))+
# For manuscript plot:
theme(
legend.position = "top",
text = element_text(size = 15.5),
axis.text.x = element_text(size = 14., margin = margin(b = 10)),
axis.text.y = element_text(size = 14.5, margin = margin(l = 10)),
legend.text = element_text(size = 17.5)
)
plot_sim_absbias_widloc_comb
ggsave(here("plots","sim_main", "sim_absbias_widloc_comb.pdf"),
plot_sim_absbias_widloc_comb, width = 11, height = 7)
```
### MSE for Location and Width
#### Combined MSE
Show both ITM and simple means/medians in a comparison:
```{r}
plot_sim_mse <- sim_res_itm_long |>
filter(pm == "mse") |>
filter(summary == "mean") |>
filter(measure %in% c("Trinterval", "simplemeaninterval", "simplemedinterval")) |>
mutate(measure = case_when(
measure == "Trinterval" ~ "ITM",
measure == "simplemeaninterval" ~ "Simple Means",
measure == "simplemedinterval" ~ "Simple Medians"
)) |>
mutate(n_items = paste0(n_items, " Items")) |>
# order n_items correctly
mutate(n_items = factor(n_items, levels = c("5 Items",
"10 Items",
"20 Items",
"40 Items"))) |>
ggplot(aes(x = n_respondents,
y = mean,
color = measure,
group = measure)) +
# add vertical line between different sample sizes
geom_vline(colour = "#F3F4F5", xintercept = seq(1.5, 4, 1))+
geom_point(position = position_dodge(0.7),
size = 2.5) +
geom_errorbar(aes(ymin = mean - 1*mcse,
ymax = mean + 1*mcse),
width = .8,
position = position_dodge(0.7),
show.legend = FALSE)+
ggh4x::facet_wrap2(n_items ~ .,
axes = "all",
nrow = 1) +
scale_y_continuous(limits = c(0, .49), expand = c(0,0)) +
ggokabeito::scale_color_okabe_ito(order = c(5, 1, 3))+
labs(x = "Number of Respondents",
y = "MSE",
color = "") +
theme_icm()+
theme(legend.position = "top",
text = element_text(size = 28))
ggsave(here("plots","sim_main", "sim_mse.pdf"),
plot_sim_mse, width = 7, height = 3.75)
plot_sim_mse
```
#### Separate MSE
Separated by location and width and standardize values with the true variance:
```{r}
plot_sim_mse_widloc <- sim_res_itm_long |>
filter(pm == "mse") |>
filter(summary == "mean") |>
filter(measure %in% c("Trloc", "Trwid")) |>
mutate(measure = case_when(measure == "Trloc" ~ "Location", measure == "Trwid" ~ "Width")) |>
mutate(n_items = paste0(n_items, " Items")) |>
# standardize with true variability
mutate(
mean = if_else(
measure == "Location",
mean / sigma_Tr_loc^2,
mean / sigma_Tr_wid^2
),
mcse = if_else(
measure == "Location",
mcse / sigma_Tr_loc^2,
mcse / sigma_Tr_wid^2
)
) |>
# order n_items correctly
mutate(n_items = factor(n_items, levels = c("5 Items", "10 Items", "20 Items", "40 Items"))) |>
ggplot(aes(
x = n_respondents,
y = mean,
color = measure,
group = measure
)) +
# add vertical line between different sample sizes
geom_vline(colour = "#F3F4F5", xintercept = seq(1.5, 4, 1)) +
geom_point(position = position_dodge(0.7), size = 2.5) +
geom_errorbar(
aes(ymin = mean - 1 * mcse, ymax = mean + 1 * mcse),
width = .8,
position = position_dodge(0.7),
show.legend = FALSE
) +
ggh4x::facet_wrap2(n_items ~ ., axes = "all", nrow = 1) +
scale_y_continuous(limits = c(0, 0.151), expand = c(0, 0)) +
ggokabeito::scale_color_okabe_ito(order = c(2, 7)) +
labs(x = "Number of Respondents", y = "MSE", color = "") +
theme_icm() +
theme(legend.position = "top", text = element_text(size = 28))
ggsave(here("plots","sim_main", "sim_mse_widloc.pdf"),
plot_sim_mse_widloc, width = 7, height = 3.75)
plot_sim_mse_widloc
```
Alternatively, combine plots for location and width and standardize values with the true variance:
```{r}
plot_sim_mse_widloc_comb <- sim_res_itm_long |>
filter(pm == "mse") |>
filter(summary == "mean") |>
filter(measure %in% c("Trloc", "Trwid", "simplemeanloc", "simplemeanwid", "simplemedloc", "simplemedwid")) |>
# split based on model
mutate(model = case_when(
grepl("^Tr", measure) ~ "ITM",
grepl("^simplemean", measure) ~ "Simple Means",
grepl("^simplemed", measure) ~ "Simple Medians"
),
measure = sub("^Tr|simplemean|simplemed", "", measure)) |>
mutate(measure = case_when(
measure == "loc" ~ "Location",
measure == "wid" ~ "Width")) |>
mutate(n_items = paste0(n_items, " Items")) |>
# standardize with true variability
mutate(mean = if_else(measure == "Location",
mean / sigma_Tr_loc^2,
mean / sigma_Tr_wid^2
),
mcse = if_else(measure == "Location",
mcse / sigma_Tr_loc^2,
mcse / sigma_Tr_wid^2
)) |>
# order n_items correctly
mutate(n_items = factor(n_items,
levels = c("5 Items", "10 Items", "20 Items", "40 Items"))) |>
ggplot(aes(
x = n_respondents,
y = mean,
color = model,
group = model
)) +
# add vertical line between different sample sizes
geom_vline(colour = "#F3F4F5", xintercept = seq(1.5, 4, 1)) +
geom_point(position = position_dodge(0.7), size = 2.5) +
geom_errorbar(
aes(ymin = mean - 1 * mcse, ymax = mean + 1 * mcse),
width = .8,
position = position_dodge(0.7),
show.legend = FALSE
) +
ggh4x::facet_grid2(measure ~ n_items,
axes = "all") +
scale_y_continuous(limits = c(0, 0.21), expand = c(0, 0)) +
ggokabeito::scale_color_okabe_ito(order = c(5, 1, 3)) +
labs(x = "Number of Respondents", y = "Standardized MSE", color = "") +
theme_icm() +
# theme(legend.position = "top", text = element_text(size = 14))+
# For manuscript plot:
theme(legend.position = "top",
text = element_text(size = 14),
axis.text.x = element_text(size = 13.5, margin = margin(b = 10)),
axis.text.y = element_text(size = 13.5, margin = margin(l = 10)),
legend.text = element_text(size = 15.5))
plot_sim_mse_widloc_comb
ggsave(here("plots","sim_main", "sim_mse_widloc_comb.pdf"),
plot_sim_mse_widloc_comb, width = 11, height = 7)
```
### Scatterplot of Location and Width Bias
Apply a function that retrieves location and width estimates to the results:
```{r}
if(exists(here(
"sim_results",
"sim_results_itm_2025-04-28_03-30-27",
"locwid_bias.rds"
))) {
# read the data
locwid_bias <- readRDS(here(
"sim_results",
"sim_results_itm_2025-04-28_03-30-27",
"locwid_bias.rds"
))
} else{
locwid_bias <- prep_locwid(path_full_results)
# save the data
saveRDS(
locwid_bias,
here(
"sim_results",
"sim_results_itm_2025-04-28_03-30-27",
"locwid_bias.rds"
)
)
}
```
Then create a scatter plot to investigate compensatory behavior, in other words, if the bias of the location is higher when the bias of the width is lower and vice versa. Overall, there does not seem to be strong evidence for compensatory behavior.
```{r}
# standardize the values
locwid_bias <- locwid_bias |>
dplyr::mutate(loc_bias_std = loc_bias / sigma_Tr_loc,
wid_bias_std = wid_bias / sigma_Tr_wid)
# if fonts aren't working properly with ggExtra
# extrafont::font_import()
# # add google font
sysfonts::font_add_google("News Cycle", "news")
# use showtext
showtext::showtext_auto()
# windows()
# extrafont::font_import(pattern = "NewsCycle-Regular.ttf")
# extrafont::loadfonts()
scatter_tmp <- locwid_bias |>
ggplot(aes(x = loc_bias_std, y = wid_bias_std)) +
geom_point(alpha = 0.5) +
geom_abline(intercept = 0,
slope = 1,
linetype = "dashed") +
labs(x = "Location Bias (standardized)", y = "Width Bias (standardized)") +
geom_smooth(method = "loess", se = TRUE) +
scale_x_continuous(limits = c(0, .85), expand = c(0, 0)) +
scale_y_continuous(limits = c(0, .85), expand = c(0, 0)) +
# manually set theme here due to font conflicts
ggplot2::theme_minimal(base_family = "news") +
ggplot2::theme(
# remove minor grid
panel.grid.minor = ggplot2::element_blank(),
# Title and Axis Texts
plot.title = ggplot2::element_text(
face = "plain",
size = ggplot2::rel(1.2),
hjust = 0.5
),
plot.subtitle = ggplot2::element_text(size = ggplot2::rel(1.1), hjust = 0.5),
axis.text.x = ggplot2::element_text(face = "plain", size = ggplot2::rel(1.05)),
axis.text.y = ggplot2::element_text(face = "plain", size = ggplot2::rel(1.05)),
axis.title.x = ggplot2::element_text(face = "plain", size = ggplot2::rel(1.3)),
axis.title.y = ggplot2::element_text(face = "plain", size = ggplot2::rel(1.3)),
axis.line = element_line(colour = "#6d6d6e"),
# Faceting
strip.text = ggplot2::element_text(
face = "plain",
size = ggplot2::rel(1.1),
hjust = 0.5
),
strip.text.x.top = ggplot2::element_text(
face = "plain",
size = ggplot2::rel(1.2),
hjust = 0.5
),
# strip.text.y = element_blank(),
strip.background = ggplot2::element_rect(fill = NA, color = NA),
# Grid
panel.grid = ggplot2::element_line(colour = "#F3F4F5"),
# Legend
legend.title = ggplot2::element_text(face = "plain"),
legend.position = "top",
legend.justification = 1,
# Panel/Facets
panel.spacing.x = ggplot2::unit(1.6, "lines"),
panel.spacing.y = ggplot2::unit(1.6, "lines"),
# Remove vertical grid lines
panel.grid.major.x = ggplot2::element_blank()
) +
theme(text = element_text(size = 28))
# show marginal distributions as histograms
# doesn't work with custom font
# scatter_locwid <- ggExtra::ggMarginal(scatter_tmp, type = "histogram", bins = 50)+
# theme_minimal(base_family = "news")
ggsave(here("plots","sim_main", "scatter_locwid.pdf"),
scatter_tmp, width = 7, height = 5)
scatter_tmp
```
### Visualization of Raw Results
Again, apply a function that retrieves location and width estimates to the results:
```{r}
if(exists(here(
"sim_results",
"sim_results_itm_2025-04-28_03-30-27",
"raw_absbias.rds"
))) {
# read the data
path_full_results <- here("sim_results",
"sim_results_itm_2025-04-28_03-30-27",
"raw_absbias.rds")
} else{
# function to obtain the raw results
raw_absbias <-
prep_locwid(path_full_results,
simplemeans = TRUE,
simplemedians = TRUE)
# save the data
saveRDS(
raw_absbias,
here(
"sim_results",
"sim_results_itm_2025-04-28_03-30-27",
"raw_absbias.rds"
)
)
}
```
Or just read in the results:
```{r}
raw_absbias <- readRDS(here("sim_results",
"sim_results_itm_2025-04-28_03-30-27",
"raw_absbias.rds"))
```
Instead of only showing aggregate performance measures, we can also visualize the raw repetition-wise performance measures to investigate the variability and potential skewness of the performance.
```{r}
# attach them to the condition-wise results
raw_ests_boxplot <- sim_res_itm |>
mutate(iteration = row_number()) |>
dplyr::select(n_respondents, n_items, iteration) |>
right_join(raw_absbias, by = "iteration") |>
pivot_longer(cols = c(loc_bias, wid_bias, smloc_bias, smwid_bias, smmloc_bias, smmwid_bias)) |>
mutate(n_respondents = factor(n_respondents)) |>
mutate(name = gsub("_bias", "", name)) |>
mutate(
name = case_match(
name,
"loc" ~ "ITM Location",
"wid" ~ "ITM Width",
"smloc" ~ "Simplemeans Location",
"smwid" ~ "Simplemeans Width",
"smmloc" ~ "Simplemedians Location",
"smmwid" ~ "Simplemedians Width",
.default = name
)
) |>
separate_wider_delim(name, delim = " ", names = c("model", "measure")) |>
mutate(model = gsub("Simplemeans", "Simple Means", model)) |>
mutate(model = gsub("Simplemedians", "Simple Medians", model)) |>
mutate(n_items = paste0(n_items, " Items")) |>
# order n_items correctly
mutate(n_items = factor(n_items, levels = c("5 Items", "10 Items", "20 Items", "40 Items"))) |>
mutate(n_respondents = as_factor(n_respondents)) |>
ggplot(aes(
x = n_respondents,
y = value,
color = model,
group = model,
fill = model
)) +
# add vertical line between different sample sizes
geom_vline(colour = "#F3F4F5", xintercept = seq(1.5, 4, 1)) +
# ggdist::stat_halfeye(
# position = "dodge",
# scale = 5
# ) +
geom_boxplot(
alpha = .2,
lwd = 0.6,
fatten = 0.9,
inherit.aes = FALSE,
aes(
x = n_respondents,
y = value,
fill = model,
color = model
)
) +
# geom_point(
# size = 1,
# alpha = 0.15,
# position = ggplot2::position_dodge(width = .6)
# )+
# geom_jitter(size=0.4, alpha=0.2) +
ggh4x::facet_grid2(measure ~ n_items, axes = "all", scales = "free_y") +
# scale_y_continuous(limits = c(0, 1), expand = c(0, 0)) +
ggokabeito::scale_color_okabe_ito(order = c(5, 1, 3)) +
ggokabeito::scale_fill_okabe_ito(order = c(5, 1, 3)) +
labs(x = "Number of Respondents",
y = "Absolute Bias",
color = "",
fill = "") +
theme_icm() +
# theme(legend.position = "top", text = element_text(size = 18))+
# For manuscript plot:
theme(
legend.position = "top",
text = element_text(size = 14),
axis.text.x = element_text(size = 13.5, margin = margin(b = 10)),
axis.text.y = element_text(size = 13.5, margin = margin(l = 10)),
legend.text = element_text(size = 15.5)
)
raw_ests_boxplot
ggsave(here("plots","sim_main", "raw_ests_boxplot.pdf"),
raw_ests_boxplot, width = 15, height = 7)
```
### Summary Table
To show the numerical results of the main outcome measures, we create a summary table below. It shows unstandardized estimates.
```{r}
sim_res_itm |>
select(
n_items,
n_respondents,
Tr_interval_mean_fn_abs_bias_mean,
Tr_interval_mean_fn_abs_bias_mcse,
Tr_interval_mean_fn_mse_mean,
Tr_interval_mean_fn_mse_mcse,
simplemean_interval_mean_fn_abs_bias_mcse,
simplemean_interval_mean_fn_abs_bias_mean,
simplemean_interval_mean_fn_mse_mean,
simplemean_interval_mean_fn_mse_mcse,
simplemed_interval_mean_fn_abs_bias_mcse,
simplemed_interval_mean_fn_abs_bias_mean,
simplemed_interval_mean_fn_mse_mean,
simplemed_interval_mean_fn_mse_mcse
) |>
pivot_longer(cols = !c(n_items, n_respondents)) |>
# split based on last underscore
separate(
name,
into = c("name", "suffix"),
sep = "_(?=[^_]+$)",
remove = FALSE
) |>
pivot_wider(names_from = suffix, values_from = value) |>
# again split based on last underscore
separate(name,
into = c("name", "pm"),
sep = "_(?=[^_]+$)",
remove = FALSE) |>
# remove "mean_fn" from name
mutate(name = str_remove(name, "_mean_fn_abs")) |>
mutate(name = str_remove(name, "_mean_fn")) |>
dplyr::rename(
"Items" = "n_items",
"Respondents" = "n_respondents",
"Model" = "name",
"Measure" = "pm",
"Mean" = "mean",
"MCSE" = "mcse"
) |>
mutate(
Model = case_when(
Model == "Tr_interval" ~ "ITM",
Model == "simplemean_interval" ~ "Simple Means",
Model == "simplemed_interval" ~ "Simple Medians"
)
) |>
mutate(Measure = case_when(Measure == "bias" ~ "Bias", Measure == "mse" ~ "MSE")) |>
arrange(Items, Respondents) |>
mutate(across(c(Mean, MCSE), ~ round(., 4))) |>
knitr::kable(align = c("c", "c", "l", "l", "c", "c"))
```
## Absolute Bias and MSE for Other Parameters
The performance measures of all other parameters can be accessed in the tabs below.
Create a plotting function:
```{r}
plot_bias <- function(data, pm, param, lab_y_axis = "Absolute Bias", ...) {
# obtain y-axis scaling
y_max <- data |>
filter(pm == {{pm}}) |>
filter(summary == "mean") |>
filter(measure == {{param}}) |>
mutate(max_val = mean + 1 * mcse) |>
pull(max_val)
y_max <- max(y_max) * 1.1
data |>
dplyr::filter(pm == {{pm}}) |>
dplyr::filter(summary == "mean") |>
dplyr::filter(measure == {{param}}) |>
dplyr::mutate(
measure = case_when(
measure == "Trloc" ~ "Consensus Location",
measure == "Trwid" ~ "Consensus Width",
measure == "Trinterval" ~ "Interval",
measure == "aloc" ~ "a Location",
measure == "bloc" ~ "b Location",
measure == "bwid" ~ "b Width",
measure == "Eloc" ~ "E Location",
measure == "Ewid" ~ "E Width",
measure == "lambdaloc" ~ "Lambda Location",
measure == "lambdawid" ~ "Lambda Width",
measure == "omegacor" ~ "Omega Correlation",
measure == "simplmeanloc" ~ "Simplemeans Location",
measure == "simplmeanwid" ~ "Simplemeans Width",
measure == "simplemeaninterval" ~ "Simplemeans Interval"
)
) |>
dplyr::mutate(n_items = paste0(n_items, " Items")) |>
# order n_items correctly
dplyr::mutate(n_items = factor(n_items, levels = c(
"5 Items", "10 Items", "20 Items", "40 Items"
))) |>
ggplot(aes(x = n_respondents, y = mean)) +
# add vertical line between different sample sizes
geom_vline(colour = "#F3F4F5", xintercept = seq(1.5, 4, 1)) +
geom_point(size = 2.5) +
geom_errorbar(
aes(ymin = mean - 1 * mcse, ymax = mean + 1 * mcse),
width = .8,
show.legend = FALSE
) +
ggh4x::facet_wrap2(n_items ~ ., axes = "all", nrow = 1) +
scale_y_continuous(expand = c(0, 0), limits = c(0, y_max)) +
# ggokabeito::scale_color_okabe_ito(order = c(2, 7))+
labs(x = "Number of Respondents", y = lab_y_axis) +
theme_icm() +
theme(legend.position = "top", text = element_text(size = 28))
}
```
### Absolute Bias
Additionally, we can plot the standardized absolute biases and MSEs on the log scale.
First, prep the data to select parameters on the log-scale:
```{r}
sim_res_log <- sim_res_itm |>
select(
n_respondents,
n_items,
contains("log"),
contains("omega"),
contains("Tr_loc"),
contains("Tr_wid"),
contains("b_loc"),
contains("b_wid")
) |>
rename_all(~ str_remove(., "_fn")) |>
pivot_longer(
cols = !c(n_respondents, n_items),
names_to = "measure",
values_to = "value"
) |>
# remove the "log", as everything is on a log scale (besides omegacor)
mutate(measure = str_remove(measure, "log_")) |>
# rename for easier separation
mutate(measure = str_replace(measure, "abs_bias", "absbias")) |>
# add string so that each measure has same number of underscores
mutate(measure = str_replace(measure, "omega", "omega_cor")) |>
# remove only the first underscore
mutate(measure = sub("_", "", measure, fixed = TRUE)) |>
separate_wider_delim(measure,
names = c("measure", "summary", "pm", "param"),
delim = "_") |>
group_by(n_respondents, n_items, measure, summary, pm) |>
pivot_wider(names_from = "param", values_from = "value") |>
ungroup() |>
mutate(n_respondents = factor(n_respondents))
```
::: {.panel-tabset}
```{r}
#| results: asis
# relevant parameters
rel_params <- c("aloc", "bloc", "bwid",
"Eloc", "Ewid", "lambdaloc", "lambdawid", "omegacor")
# set names to provide proper names of parameters
rel_param_names <- c(
aloc = "a Location (log-scale)",
bloc = "b Location",
bwid = "b Width (log-scale)",
Eloc = "E Location (log-scale)",
Ewid = "E Width (log-scale)",
lambdaloc = "Lambda Location (log-scale)",
lambdawid = "Lambda Width (log-scale)",
omegacor = "Omega Correlation"
)
# create tabs
for (param in rel_params) {
cat("###", param, "\n\n")
cat("\n \n")
print(
plot_bias(sim_res_log, pm = "absbias", param = param) +
labs(
title = rel_param_names[[param]]
)
)
cat(' \n \n')
}
```
:::
### MSE
::: {.panel-tabset}
```{r}
#| results: asis
# relevant parameters
rel_params <- c("aloc", "bloc", "bwid",
"Eloc", "Ewid", "lambdaloc", "lambdawid", "omegacor")
# set names to provide proper names of parameters
rel_param_names <- c(
aloc = "a Location (log-scale)",
bloc = "b Location",
bwid = "b Width (log-scale)",
Eloc = "E Location (log-scale)",
Ewid = "E Width (log-scale)",
lambdaloc = "Lambda Location (log-scale)",
lambdawid = "Lambda Width (log-scale)",
omegacor = "Omega Correlation"
)
# create tabs
for (param in rel_params) {
cat("###", param, "\n\n")
cat("\n \n")
print(
plot_bias(sim_res_log, pm = "mse", param = param, lab_y_axis = "MSE") +
labs(
title = rel_param_names[[param]]
)
)
cat(' \n \n')
}
```
:::
## Correlation of Estimates with True Parameters
Define the plotting function for correlations
```{r}
plot_corr <- function(data, pm, param, ...) {
# obtain y-axis scaling
y_max <- data |>
filter(pm == {{pm}}) |>
filter(summary == "mean") |>
filter(measure == {{param}}) |>
mutate(max_val = mean + 1 * mcse) |>
pull(max_val)
y_max <- max(y_max) + .005
y_min <- data |>
filter(pm == {{pm}}) |>
filter(summary == "mean") |>
filter(measure == {{param}}) |>
mutate(min_val = mean - 1 * mcse) |>
pull(min_val)
y_min <- min(y_min) - .005
# plot
data |>
dplyr::filter(pm == {{pm}}) |>
dplyr::filter(summary == "mean") |>
dplyr::filter(measure == {{param}}) |>
dplyr::mutate(
measure = case_when(
measure == "Trloc" ~ "Consensus Location",
measure == "Trwid" ~ "Consensus Width",
measure == "aloc" ~ "a Location (log-scale)",
measure == "bloc" ~ "b Location",
measure == "bwid" ~ "b Width",
measure == "Eloc" ~ "E Location (log-scale)",
measure == "Ewid" ~ "E Width (log-scale)",
measure == "lambdaloc" ~ "Lambda Location (log-scale)",
measure == "lambdawid" ~ "Lambda Width (log-scale)",
measure == "omegacor" ~ "Omega Correlation"
)
) |>
dplyr::mutate(n_items = paste0(n_items, " Items")) |>
# order n_items correctly
dplyr::mutate(n_items = factor(n_items, levels = c(
"5 Items", "10 Items", "20 Items", "40 Items"
))) |>
ggplot(aes(x = n_respondents, y = mean)) +
# add vertical line between different sample sizes
geom_vline(colour = "#F3F4F5", xintercept = seq(1.5, 4, 1)) +
geom_hline(yintercept = 1,
linetype = "dashed",
color = "grey50") +
geom_point(size = 2.5) +
geom_errorbar(
aes(ymin = mean - 1 * mcse, ymax = mean + 1 * mcse),
width = .8,
show.legend = FALSE
) +
ggh4x::facet_wrap2(n_items ~ ., axes = "all", nrow = 1) +
scale_y_continuous(expand = expansion(), limits = c(y_min, y_max)) +
# ggokabeito::scale_color_okabe_ito(order = c(2, 7))+
labs(x = "Number of Respondents", y = "Correlation") +
theme_icm() +
theme(legend.position = "top", text = element_text(size = 28))
}
```
::: {.panel-tabset}
```{r}
#| results: asis
# relevant parameters
rel_params <- c(
"Trloc",
"Trwid",
"aloc",
"bloc",
"bwid",
"Eloc",
"Ewid",
"lambdaloc",
"lambdawid",
"omegacor"
)
# set names to provide proper names of parameters
rel_param_names <- c(
Trloc = "Consensus Location",
Trwid = "Consensus Width",
aloc = "a Location (log-scale)",
bloc = "b Location",
bwid = "b Width",
Eloc = "E Location (log-scale)",
Ewid = "E Width (log-scale)",
lambdaloc = "Lambda Location (log-scale)",
lambdawid = "Lambda Width (log-scale)",
omegacor = "Omega Correlation"
)
# create tabs
for (param in rel_params) {
cat("###", param, "\n\n")
cat("\n \n")
print(
plot_corr(sim_res_log, pm = "corr", param = param) +
labs(
title = rel_param_names[[param]]
)
)
cat(' \n \n')
}
```
:::
### All Parameters in One Plot
```{r}
# min for y-axis
y_min <-
sim_res_log |>
dplyr::filter(pm == "corr") |>
pull(mean) |>
min()
# data
plot_corr_all_pars <- sim_res_log |>
dplyr::filter(pm == "corr") |>
dplyr::filter(summary == "mean") |>
dplyr::mutate(
measure = case_when(
measure == "Trloc" ~ "T^loc",
measure == "Trwid" ~ "T^wid",
measure == "aloc" ~ "a^loc~'(log)'",
measure == "bloc" ~ "b^loc",
measure == "bwid" ~ "b^wid",
measure == "Eloc" ~ "E^loc~'(log)'",
measure == "Ewid" ~ "E^wid~'(log)'",
measure == "lambdaloc" ~ "lambda^loc~'(log)'",
measure == "lambdawid" ~ "lambda^wid~'(log)'",
measure == "omegacor" ~ "omega"
)
) |>
dplyr::mutate(n_items = paste0(n_items, " Items")) |>
# order n_items correctly
dplyr::mutate(n_items = factor(n_items, levels = c("5 Items", "10 Items", "20 Items", "40 Items"))) |>
# plot
ggplot(aes(
x = n_respondents,
y = mean,
color = n_items,
group = n_items
)) +
# add vertical line between different sample sizes
geom_vline(colour = "#F3F4F5", xintercept = seq(1.5, 4, 1)) +
geom_point(position = position_dodge(0.7), size = 2.5) +
geom_errorbar(
aes(ymin = mean - 1 * mcse, ymax = mean + 1 * mcse),
width = .8,
position = position_dodge(0.7),
show.legend = FALSE
) +
scale_y_continuous(
breaks = seq(0.5, 1, .1),
labels = c(".5", ".6", ".7", ".8", ".9", "1"),
limits = c(y_min - 0.02, 1.02),
expand = expansion()
) +
facet_wrap2(
vars(measure),
ncol = 5,
axes = "all",
scales = "free",
labeller = labeller(measure = label_parsed)
) +
ggokabeito::scale_color_okabe_ito(order = c(7, 5, 1, 3)) +
labs(x = "Number of Respondents", y = "Correlation", color = NULL) +
theme_icm() +
# theme(legend.position = "top", text = element_text(size = 14))+
# For manuscript plot:
theme(
legend.position = "top",
text = element_text(size = 14),
axis.text.x = element_text(size = 14, margin = margin(b = 10, t = 10)),
axis.title.x = element_text(margin = margin(t = 12)),
axis.text.y = element_text(size = 14),
legend.text = element_text(size = 16),
strip.background = element_blank(),
strip.text = element_text(margin = margin(0, 0, 0, 0))
)
plot_corr_all_pars
ggsave(
filename = here("plots", "sim_main", "sim_corr.pdf"),
plot = plot_corr_all_pars,
width = 11,
height = 8,
scale = 1
)
```
## All Measures for All Parameters in One Summary Table
```{r}
sim_res_log |>
dplyr::filter(measure %in% rel_params,
summary == "mean",
pm %in% c("absbias", "mse", "corr")) |>
# round mean and mcse
dplyr::mutate(mean = round(mean, 3), mcse = round(mcse, 3)) |>
# pivot_wider of "pm"
pivot_wider(names_from = pm, values_from = c(mean, mcse, sd)) |>
# reorder columns
select(
n_respondents,
n_items,
measure,
mean_absbias,
mcse_absbias,
mean_mse,
mcse_mse,
mean_corr,
mcse_corr
) |>
# replace values in measure with proper names
mutate(
measure = case_when(
measure == "Trloc" ~ "Consensus Location",
measure == "Trwid" ~ "Consensus Width",
measure == "aloc" ~ "a Location (log-scale)",
measure == "bloc" ~ "b Location",
measure == "bwid" ~ "b Width",
measure == "Eloc" ~ "E Location (log-scale)",
measure == "Ewid" ~ "E Width (log-scale)",
measure == "lambdaloc" ~ "Lambda Location (log-scale)",
measure == "lambdawid" ~ "Lambda Width (log-scale)",
measure == "omegacor" ~ "Omega Correlation",
TRUE ~ measure
)
) |>
rename(
"Items" = n_items,
"Respondents" = n_respondents,
Parameter = measure,
"Mean Absolute Bias" = mean_absbias,
"MCSE Absolute Bias" = mcse_absbias,
"Mean MSE" = mean_mse,
"MCSE MSE" = mcse_mse,
"Mean Correlation" = mean_corr,
"MCSE Correlation" = mcse_corr
) |>
mutate(
Items = factor(Items),
Respondents = factor(Respondents),
Parameter = factor(Parameter)
) |>
arrange(Items, Respondents) |>
#knitr::kable(align = c("c", "c", "l", "c", "c", "c", "c"))
DT::datatable(
filter = "top",
rownames = FALSE,
options = list(
pageLength = 50,
lengthMenu = c(50, 100),
dom = "t",
columnDefs = list(list(
className = 'dt-center',
targets = c(
"Items",
"Respondents",
"Mean Absolute Bias",
"MCSE Absolute Bias",
"Mean MSE",
"MCSE MSE",
"Mean Correlation",
"MCSE Correlation"
)
)),
autoWidth = TRUE,
scrollX = TRUE,
scrollY = "720px"
)
)
```
***
# Session Info
```{r}
pander::pander(sessionInfo())
```
